#ifndef LUAEIGEN_VECTOR_HPP
#define LUAEIGEN_VECTOR_HPP



#include <luaNumeric/config.h>
#include <vector>
#include <string>
#include <Eigen/Core>
#include <Eigen/Dense>

#include <boost/bind.hpp>
#include <boost/math/special_functions/sin_pi.hpp>
#include <algorithm>
#include <numeric>
#include <complex>


namespace Vector{
	std::string serialize( const std::vector<double>& vec );
	std::string serialize( const std::vector<std::complex<double>>& vec );

	template<class T>
	T min( const std::vector<T>& vec )
	{
		return *std::min_element(vec.begin(),vec.end());
	}

	template<class T>
	void min_and_index( const std::vector<T>& vec, T& val, int& idx )
	{
		auto ele = std::min_element(vec.begin(),vec.end());
		val = *ele;
		idx = std::distance(vec.begin(),ele);
	}

	template<class T>
	void min_and_index_with_range( const std::vector<T>& vec, int min, int max, T& val, int& idx )
	{
		if( min > max ) std::swap(min,max);
		min = std::min( std::max( 0, min), (int)vec.size() );
		max = std::max( std::min( (int)vec.size(), max ), 0 );
		auto it1 = vec.begin()+min;
		auto it2 = vec.begin()+max;
		auto ele = std::min_element(it1,it2);
		val = *ele;
		idx = std::distance(vec.begin(),ele);
	}

	template<class T>
	void max_and_index( const std::vector<T>& vec, T& val, int& idx )
	{
		auto ele = std::max_element(vec.begin(),vec.end());
		val = *ele;
		idx = std::distance(vec.begin(),ele);
	}

	template<class T>
	void max_and_index_with_range( const std::vector<T>& vec, int min, int max, T& val, int& idx )
	{
		if( min > max ) std::swap(min,max);
		min = std::min( std::max( 0, min), (int)vec.size() );
		max = std::max( std::min( (int)vec.size(), max ), 0 );
		auto it1 = vec.begin()+min;
		auto it2 = vec.begin()+max;
		auto ele = std::max_element(it1,it2);
		val = *ele;
		idx = std::distance(vec.begin(),ele);
	}

	template<class T>
	T interpolation(const std::vector<double>& xarr, const std::vector<T>& yarr,double x)
	{
		auto pos = std::find_if(xarr.begin(),xarr.end(),boost::bind(std::greater<double>(),_1,x));
		int npos = std::distance(xarr.begin(),pos);
		int nlow = std::max(npos-1,0);
		int nhigh = std::min((int)xarr.size()-1,npos);
		if( npos == 0) return yarr[0];
		if( npos == xarr.size() ) return yarr[xarr.size()-1];
		T res = yarr[npos-1] + (x-xarr[npos-1])/(xarr[npos]-xarr[npos-1])*(yarr[npos]-yarr[npos-1]);
		return res;
	}


	template<class T>
	std::vector<T> interpolation(const std::vector<double>& xarr, const std::vector<T>& yarr,const std::vector<double>& x)
	{
		std::vector<T> y(x.size());
		for( int ii=0; ii < y.size(); ii++){
			y[ii] = interpolation(xarr,yarr,x[ii]);
		}
		return y;
	}
	template<class T>
	std::vector<double> abs( const std::vector<T>& vec );
	template<class T>
	std::vector<double> phase( const std::vector<T>& vec );

	template<class T>
	std::vector<T> interpolation_complex(const std::vector<double>& xarr,
						                 const std::vector<T>& yarr,const std::vector<double>& x)
	{
		auto amp = Vector::abs(yarr);
		auto phs = Vector::phase(yarr);
		auto amp1 = interpolation(xarr,amp,x);
		auto phs1 = interpolation(xarr,phs,x);

		std::vector<T> y(x.size());
		for( int ii=0; ii < y.size(); ii++){
			y[ii] = std::polar(amp1[ii], phs1[ii]);
		}
		return y;
	}

	template<class T>
	void sort( std::vector<T>& vec )
	{
		std::sort(vec.begin(),vec.end());
	}

	template<class CONT>
	bool grids_index_impl(
					 //std::vector< std::vector<int> >& sheet,
			 	 	 CONT& sheet,
			         int col_idx,			  //当前操作列
					 int nrows,               //当前列的数量
					 int ncols                //变量数量
					 )
	{
		//
		if(nrows<1){
			throw("index out of range at grids operation");
			return false;
		}

		if( !sheet.size() ){
			for(int ii=0; ii < nrows; ii++ ){
				std::vector<int> data(ncols);
				data[col_idx] = ii;
				sheet.push_back(data);
			}
			return true;
		}

		int nold = sheet.size();

		for( int jj=0; jj < nold; jj++ ){
			sheet[jj][col_idx] = 0;
		}

		for( int ii=1; ii < nrows; ii ++){
			for( int jj=0; jj < nold; jj++ ){
				std::vector<int> data = sheet[jj];
				data[col_idx] = ii;
				sheet.push_back(data);
			}
		}
		return true;
	}

	template<class T>
	bool grids_index( const std::vector< std::vector< T > >& vec, std::vector< std::vector<int> >& idx )
	{
		if(!vec.size()) return false;

		int ncols = vec.size();
		int nlists = 1;
		for(auto& v:vec){
			nlists *= v.size();
		}

		idx.clear();
		idx.reserve(nlists);


		for(int ii=0; ii < vec.size(); ii++){
			if(! grids_index_impl(idx,ii,vec[ii].size(),ncols) )
				return false;
		}
		return true;
	}

	template<class T>
	std::tuple< std::vector< T >,std::vector< T > > grids( const std::vector< T >& x,
			          const std::vector< T >& y )
	{
		int xnum = x.size();
		int ynum = y.size();
		int nn = xnum*ynum;
		std::vector<T> xx(nn),yy(nn);
		for( int ii=0; ii < xnum; ii ++)
		{
			for( int jj=0; jj < ynum; jj ++)
			{
				xx[ii*ynum+jj] = x[ii];
				yy[ii*ynum+jj] = y[jj];
			}
		}
		return std::make_tuple(std::move(xx),std::move(yy));
	}


	template<class T>
	T max( const std::vector<T>& vec )
	{
		return *std::max_element(vec.begin(),vec.end());
	}
	template<class T>
	T sum( const std::vector<T>& vec )
	{
		T res(0);
		return std::accumulate(vec.begin(),vec.end(),res);
	}

	template<class T>
	std::vector<T> log10( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::log10(res[ii]);
		return res;
	}

	template<class T>
	std::vector<T> log( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++){
			res[ii] = std::log(res[ii]);
		}
		return res;
	}
	template<class T>
	std::vector<T> exp( const std::vector<T>& vec)
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::exp(res[ii]);
		return res;
	}

	template<class T>
	std::vector<T> sqrt( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::sqrt(res[ii]);
		return res;
	}

	template<class T>
	std::vector<T> sin( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::sin(res[ii]);
		return res;
	}
	template<class T>
	std::vector<T> sinc_pi( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = boost::math::sinc_pi(res[ii]);
		return res;
	}

	template<class T>
	std::vector<T> cos( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::cos(res[ii]);
		return res;
	}

	template<class T>
	std::vector<T> asin( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++){
			res[ii] = std::asin(res[ii]);
		}
		return res;
	}

	template<class T>
	std::vector<T> acos( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::acos(res[ii]);
		return res;
	}
	template<class T>
	std::vector<T> atan( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::atan(res[ii]);
		return res;
	}

	template<class T>
	std::vector<T> tan( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::tan(res[ii]);
		return res;
	}

	template<class T>
	std::vector<T> sinh( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::sinh(res[ii]);
		return res;
	}

	template<class T>
	std::vector<T> cosh( const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::cosh(res[ii]);
		return res;
	}



	template<class T>
	std::vector<T> pow( const std::vector<T>& vec, double e )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::pow(res[ii],e);
		return res;
	}

	template<class T>
	std::vector<T> pow( double e, const std::vector<T>& vec )
	{
		std::vector<T> res = vec;
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::pow(e,res[ii]);
		return res;
	}
	template<class T>
	std::vector<T> vsgn( const std::vector<T>& vec )
	{
		std::vector<T> res(vec.size(),0);
		for( int ii=0; ii < res.size(); ii++){
			if(vec[ii]>0){
				res[ii] = 1;
			}
			if(vec[ii]<0){
				res[ii] = -1;
			}
		}
		return res;
	}
	template<class T>
	T sgn( T in )
	{
		if(in>0){
			return 1;
		}
		if(in<0){
			return -1;
		}
		return 0;
	}


	template<class T>
	T dot( const std::vector<T>& vec1, const std::vector<T>& vec2 )
	{
		if( vec1.size() != vec2.size() )
			throw ( "The size of vectors dismatch!" );
		T res(0);
		for( int ii=0; ii < vec1.size(); ii++)
			res += vec1[ii]*vec2[ii];
		return res;
	}

	template<class T>
	std::vector<double> abs( const std::vector<T>& vec )
	{
		std::vector<double> res(vec.size());
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::abs(vec[ii]);
		return res;
	}

	template<class T>
	std::vector<double> arg( const std::vector<T>& vec )
	{
		std::vector<double> res(vec.size());
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::arg(vec[ii]);
		return res;
	}

	std::vector<double> arg( const std::vector<double>& vecx,const std::vector<double>& vecy );


	template<class T>
	std::vector<double> phase( const std::vector<T>& vec )
	{
		std::vector<double> res(vec.size());
		double phs0;
		for( int ii=0; ii < res.size(); ii++)
		{
			if( ii == 0 )
			{
				phs0 = std::arg( vec[ii] );
			}
			else
			{
				phs0 = phs0 + std::arg(vec[ii]/vec[ii-1]);
			}
			res[ii] = phs0;
		}
		return res;
	}

	template<class T>
	std::vector<double> real( const std::vector<T>& vec )
	{
		std::vector<double> res(vec.size());
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::real(vec[ii]);
		return res;
	}

	template<class T>
	std::vector<double> imag( const std::vector<T>& vec )
	{
		std::vector<double> res(vec.size());
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::imag(vec[ii]);
		return res;
	}
	template<class T>
	std::vector<T> conj( const std::vector<T>& vec )
	{
		std::vector<T> res(vec.size());
		for( int ii=0; ii < res.size(); ii++)
			res[ii] = std::conj(vec[ii]);
		return res;
	}

	template<class VT, class T>
	void fit_line(const VT& xx, const VT& yy, T& la, T& lb)
	{
		if( xx.size() != yy.size() ) throw ("size dismatch");
		Eigen::MatrixXd Al(xx.size(),2) ;
		Eigen::VectorXd b(yy.size());

		for(int ii=0; ii < xx.size(); ii++){
			Al(ii,0) = xx[ii];
			Al(ii,1) = 1.0;
			b(ii) = yy[ii];
		}
		Eigen::MatrixXd res = Al.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
		la = res(0);
		lb = res(1);
	}

	template<class VT>
	Eigen::MatrixXd FitterLeastSquareMethod(const VT& xx, const VT& yy, int orders)
	{
		// abnormal input verification
		if (xx.size() < 2 || yy.size() < 2 || xx.size() != yy.size() || orders < 1)
			throw "size dismatch";

		Eigen::MatrixXd Al(xx.size(),orders+1) ;
		Eigen::VectorXd b(yy.size());

		for(int ii=0; ii < xx.size(); ii++){
			for(int jj=0; jj <= orders; jj++){
				if(jj==0) {
					Al(ii,orders-jj) = 1;
					continue;
				}
				Al(ii,orders-jj) = Al(ii,orders-jj+1) * xx[ii];
			}
			b(ii) = yy[ii];
		}
		Eigen::MatrixXd res = Al.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
		return res;
	}

	template<class VT, class T>
	void fit_plane(const VT& xx, const VT& yy, const VT& zz, T& la, T& lb, T& lc, T& ld)
	{
		if (xx.size() < 2 || xx.size() != yy.size()
				|| xx.size() != zz.size() )
			throw "size dismatch";

		T xave = Vector::sum(xx)/xx.size();
		T yave = Vector::sum(yy)/yy.size();
		T zave = Vector::sum(zz)/zz.size();

		Eigen::MatrixXd Al(xx.size(),4) ;
		Eigen::VectorXd b(yy.size());

		for(int ii=0; ii < xx.size(); ii++){
			Al(ii,0) = xx[ii]-xave;
			Al(ii,1) = yy[ii]-yave;
			Al(ii,2) = zz[ii]-zave;
			Al(ii,3) = 1;
			b(ii) = 0;
		}

		Eigen::MatrixXd res = Al.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
		la = res(0);
		lb = res(1);
		lc = res(2);
		//ld = res(3);
	}

	template<class VT, class T>
	void fit_parabolic_line(const VT& xx, const VT& yy, T& la, T& lb, T& lc)
	{
		auto res = FitterLeastSquareMethod(xx,yy,2);
		la = res(0);
		lb = res(1);
		lc = res(2);
	}

	template<class T>
	std::vector<T> linear_space(const T& start, const T& end, int num){
		if( num < 1){
			return std::vector<T>();
		}else if(num==1){
			return std::vector<T>(1,start);
		}
		std::vector<T> res(num);
		T dx = (end-start)/(num-1);
		for(int ii=0; ii < num; ii++){
			res[ii] = start + ii * dx;
		}
		return std::move(res);
	}

};

#endif
